
SPECIFICATION 



TITLE 



COMPUTER-AIDED SIMULATION METHOD FOR DETERMINING THE 
ELECTROMAGNETIC FIELD OF A BODY 



BACKGROUND OF THE INVENTION 



The invention relates to a computer-aided simulation method 
.for. determining... the .electromagnetic field of a body which 
comprises a plurality of subregions and contains a plurality of 
charges and currents. 

DIN standard VDE 0870 describes the "electromagnetic 
compatibility" (EMV) of an electrical device as the "capacity of 
an electrical device to function satisfactorily in its 
electromagnetic environment without unacceptably affecting this 
environment, which is shared by other devices". 

An algorithm which belongs to the CG (conjugate gradient) 
method class and is referred to as the GMRES method, is disclosed 
by "Youcef Saas and Martin H. Schultz: GMRES: A Generalized 
Minimal Residual Algorithm for Solving Nonsymmetric Linear 
Systems. SIAM J. Sci . Stat. Comp . , Vol. 7, No. 3, pp. 856-869, 
July 1986" . 

The fast multipole method (FMM) is disclosed, for example, 
by "V. Rokhlin: Rapid Solution of Integral Equations of Classical 
Potential Theory. Journal of Computational Physics, Vol. 60, pp. 
187-207, 1985." One essential disadvantage with this formalism 
is due to the property that the multipole coeef f icients are not 
calculated explicitly, which may result in spurious distortions 
that in turn lead to errors in the multipole expansions. 
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SUMMARY OF THE INVENTION 
The object of the invention is to determine the 
electromagnetic field of a body using a computer, so that 
electromagnetic compatibility of the body can be ensured actually 
before it is made, and thus to avoid optimization cycles 
consisting of measurements and retrospective improvements, the 
above-mentioned disadvantages being avoided. 

. In general terms,, the present .invention.. j.s,...a ..comp.uter.-ai_ded 
simulation method for determining the electromagnetic field of a 
body which has a plurality of subregions and contains a plurality 
of charges and currents. In each of the plurality of subregions, 
a global multipole expansion is made which represents the effect 
of the charges and currents for distant points in the form of a 
multipole expansion, and a local multipole expansion is made, 
which represents the effect of the charges and currents at points 
inside this one of the plurality of subregions in the form of a 
multipole expansion. The electromagnetic field of the body is 
determined by superposition using the global multipole expansion 
and the local multipole expansion for the plurality of 
subregions. 

Advantageous developments of the present invention are as 
follows . 

The following steps for determining the electromagnetic 
field of the body are carried out iteratively until the error 
measure is of a predeterminable small size, 1=0 being taken as 
an initial condition: 

a) calculating the global multipole expansions with global 
multipole coefficients according to 
c9 = GI , 

c 3 being a vector made up of the global multipole 
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coefficients of the plurality of subregions, 
I being a vector which specifies a given current 
distribution, 

G being a matrix determining the global multipole 

coefficients in the relevant subregion of the plurality 
of subregions for the given current distribution I ; 

b) calculating the local multipole expansion with local 
multiple coefficients according to 

c 1 = Tcs, 

c 1 being a vector made up of the local multipole 
coefficients of the plurality of subregions, 

T being a translation matrix through which the global 
multipoles are combined into local multipoles ,- 

c) determining the electromagnetic field from 
ZI = Z'l + Lc 1 , 

Z denoting an impedance matrix, 

Z' denoting a part of the impedance matrix Z, representing 

couplings between the subregions, 
L denoting a matrix for evaluating the local multipole 

coefficients . 

In one embodiment, the subregions are of equal size. 

The size of the subregions are proportional to the distance 
from an observer region. 

The relevant subregion of the plurality of subregions is in 
each case assigned to a zone with uniform physical attribute. 

In a refinement step, the subregion of the plurality of 
subregions is split into up to eight zones. 

An element which has an impedance and is a component of the 
subregion of the body is taken into account directly in the 
matrix Z' as an impedance. 
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The electromagnetic field is determined for predeterminable 
frequencies . 

The predeterminable frequencies are determined by a minimum 
frequency and by a maximum frequency, the method being started at 
the minimum frequency and the electromagnetic field being in each 
case determined, continuing as far as the maximum frequency, with 
a predeterminable step size. 

... The .predeterminable frequencies are determined by a minium 
frequency and by a maximum frequency, the method being started at 
the maximum frequency and the electromagnetic field being in each 
case determined, continuing as far as the minimum frequency, with 
a predeterminable step size. 

The predeterminable frequencies are deterined by a minimum 
frequency and by a maximum frequency, the method being started at 
a frequency between the minimum frequency and the maximum 
frequency, and the electromagnetic field being in each case 
determined, continuing as far as the maximum frequency or as far 
as the minimum frequency, with a predeterminable step size. 

Stability at low frequencies is ensured by carrying out the 
global multipole expansions by elements instead of using basis 
functions. 

The electromagnetic compatibility of the body is determined. 

In a computer-aided simulation method, the electromagnetic 
field of a body which comprises a plurality of subregions and 
contains a plurality of charges - pole expansion and a local 
multipole expansion are made in the subregions, a global 
multipole expansion representing the effect of the charges and 
currents for distant points in the form of a multipole expansion 
and, in similar fashion, a local multipole expansion representing 
the effect of the charges and currents at points inside one of 
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the plurality of subregions in the form of a multipole expansion. 
Global multipole coefficients are calculated as parameters of the 
global multipole expansion, and local multipole coefficients are 
calculated as parameters of the local multipole expansion. 
Finally, the electromagnetic field of the body is determined from 
the global and local multipole expansions. 

Simulating the prescribed measurements needed for EMV by 
using- numerical methods substantially facilitates •optimization of 
the product components. The repeated process of measurement and 
retrospective improvement, the simulated measurement being 
possible at all points, and therefore associated costs and time 
fan be substantially reduced with the method according to the 
invention. 

The method proposed here is suitable, in particular, for 
simulating scattering problems, in which the geometrical 
structures of the models, and therefore the dimensions of the 
subregions as well, are smaller than the wavelength. Under such 
conditions, the interference effects remain limited, so that just 
a few multipole coefficients are enough to describe potentials or 
field strengths. With the stabilization method described here, 
the method is table even for low frequencies. 

A first development of the invention consists in dividing 
the body into at least one further subregion. The body therefore 
consists of subregions. 

A further development consists in carrying out the following 
steps for determining the electromagnetic field iteratively with 
an initial condition 1=0 until a truncation condition which 
gives an error measure sufficiently small size. 

A subsequent development consists in choosing subregions of 
equal size. 
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In a development, the size of the subregions is chosen to be 
proportional to the distance from an observer region. 

Yet another development consists in respectively assigning a 
subregion to a zone with uniform physical attribute. This may, 
for example, involve large surfaces with an equal electromagnetic 
property, for example equal charge distribution, or the like. 

In the context of another development, components which have 
an impedance and lie inside the body in question are taken into 
account directly as an impedance in the matrix Z' . 

The electromagnetic field may furthermore be determined for 
predeterminable frequencies . 

An additional development consists in choosing the 
predeterminable frequencies from a lowest frequency to be 
examined up to a highest frequency to be examined, continuously 
with a predeterminable step size, in order to determine the 
electromagnetic field. In this case, it is possible to start 
with the lowest frequency to be examined and continue to the 
highest frequency to be examined, or vice versa. It is also 
possible to start between the lowest and highest frequencies to 
be examined, the method for determining the electromagnetic field 
being continued to the lowest or to the highest frequency. 

An additional development consists in ensuring the stability 
of the method at low frequencies by carrying out the global 
multipole expansions by elements instead of using basis 
functions . 

One development consists in using the method to determine 
the electromagnetic compatibility of a body. 
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BRIEF DESCRIPTION OF THE DRAWINGS 
The features of the present invention which are believed to 
be novel, are set forth with particularity in the appended 
claims. The invention, together with further objects and 
advantages, may best be understood by reference to the following 
description taken in conjunction with the accompanying drawings, 
in the several Figures of which like reference numerals identify 
li-ke- elements, -and in -which: 

Fig. 1 shows a block diagram showing the steps of a method for 
determining the electromagnetic field of a body; 

Fig. 2 shows a perfectly conducting body in an external field; 

Fig. 3 shows geometrical dimensions for a line weighting 
method; 

Fig. 4 shows the effect of an n-th element pair on an m-th 
element pair; 

Fig. 5a shows a distributed impedance for elements which have 
impedance ;. 

Fig. 5b shows a concentrated impedance for elements which have 
impedance ; 

Fig. 6 shows the division of the bodywork of a car into 
subregions for the refinement stages 2, 3, and 4; 

Fig. 7 shows the division into subregions in the case of 
varying element sizes; 

Fig. 8 shows a sketch representing the spherical coordinates 
r , 6 , a ; 

Fig. 9a shows a sketch representing the geometrical dimensions 

for the global multipole expansion; 
Fig. 9b shows a sketch representing the geometrical dimensions 

for the local multipole expansion; 
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Fig. 10a shows a sketch representing a conductor loop and the 
division of the basis functions into four subregions, 

Fig. 10b shows a sketch showing the charges which are formed at 
the edges of the subregions, cannot be resolved further 
and lead to truncation errors in the expansion of 

Fig. 11 shows a sketch representing the directly and indirectly 
coupled subregions; 

-Fig. -12 •• -shows a -sketch representing .the geometrical dimensions 
for the translation; and 

Fig. 13 shows a sketch representing the directly and indirectly 



DESCRIPTION OF THE PREFERRED EMBODIMENTS 
The invention will firstly be explained in more detail with 
reference to Figures 2 to 13 . 



The Maxwell equations in differential form, as are 
sufficiently well-known to the person skilled in the art, will be 
used as a starting point : 



coupled subregions for a two- stage FMM. 



THE RETARDED POTENTIALS 



0 



(la) 



V x E = -ico/x 0 H 



(lb) 



V 



H 



= 0 



(lc) 



V x H = J + icoe 0 E 



(Id) 
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V denoting the del operator, 

H denoting the magnetic field strength, 

E denoting the electric field strength, 

p denoting the space charge density, 

e 0 denoting the permittivity of a vacuum, 

i denoting an imaginary number (V"-l) , 

to denoting the angular frequency, 

- Ho ■ denoting the permeability -of -a vacuum, 

J denoting the volume current density. 



A continuity equation follows directly from the Maxwell 

equations by combining equation (la) with divergence of equation 
(Id) : 

V • J = -lojp (2) 
If the potentials are set up as follows 

E = -icoX - Vcp (3a) 



H = 1 V x A (3b) 
Mo 



then, taking into account the Lorentz condition 

V • A + icoe 0 ^ 0 (p = 0 (4) 

this gives the two differential equations: 

(v 2 + k 2 )a = -no J (5a) 

and 
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^V 2 + k 2 j (p = - (5b) 



with the wave number 



CO 271 



^ c 0 K (6) , 

whereby 

A denoting a magnetic vector potential, 

cp denoting an electric scalar potential, 

c 0 denoting the velocity of light in a vacuum, 

A denoting the wavelength. 



Both the scalar electric potential 0 and the three Cartesian 
components of the magnetic vector potential A are described by 
the same type of differential equation which, amongst other 
things, is termed the Helmholtz equation: 

(v 2 + k 2 ) M/(r) = f (r) (7) 

ijj denoting a basis function, 

r denoting a position vector, 

f () denoting the function of an argument. 



If the solution of the Helmholtz equation is known for a 
Dirac- function excitation on the right-hand side, then the total 
solution for a given function f (r) can be determined through 
superposition. On symmetry grounds, a Dirac -function excitation 
at the origin gives rise to a spherically symmetric solution 
function , so that equation (7) in spherical coordinates can be 
reduced to the following differential equation as a function of 
r . 
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1 d o d o ? 

— — r 2 — V + kV 2 = 0 

r 2 dr dr (8) 
If the substitution 

vj> = rv|/ 

is made in equation (8) , this results in an ordinary differential 
equation with constant coefficients: 



\p" + k 2 vj> = 0 (9) 



with the general solution: 

q> = c x e- ikr + c 2 e ikr {10) 

or 



-ikr -ikr 

M/ = Cl ^ + C 2 ^— (ID 



The first term in equation (11) describes an emergent wave, 
and the second term an incident wave. Since the latter is ruled 
out on physical grounds, all that remains is to determine the 
constant C 1 . For the electric scalar potential cp, this is given 
as follows by comparing the general solution with the 
electrostatic potential (k-*0) of a charge element relative to 
pdV located at the origin: 

C? 1 = -^pdV (12) 

4TCEQ 

Similarly, the constant for the magnetic vector potential of 
a current element JdV located at the origin is given as: 

c mag = M Jdv (13) 



471 
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Through integration over all the source regions G, the 
solutions of equations (5a) and (5b) are finally given as: 



3 ; 



I: 1 



— ikR 

= ^/JJ J W^-^' (14a) 
G 

U G 

with: R = |(r — r'|| 



whereby 

R denoting the distance between an observer point r and a 

source point r' , 
G denoting a volume region, 
A, <p denoting retarded potentials. 

CALCULATION OF THE FIELD STRENGTHS 

Substituting equation (14a) in (3a) and equation (14b) in 
(3b) makes it possible to determine the electric field strength 
and the magnetic field strength: 

-ikR , -ikR 

_ xco^ rrr j(rl) e dy , _ _±_ jjj (r , )v ^ 

471 G ^ ^ G R (15a) 



-ikR 



H(r) = — fff V - x J(r') dV 

v ' 471 JJJ R • 



(15b) 
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Using spherical coordinates gives the following: 



-ikR 



= -e 



-ikR 



W 3 



R 2 ^ 



R 



(16) 



The integral formulae for determining the electric and 
magnetic field strengths for given sources are given as: 



(17a) 



(17b) 



Equation (17b) can be regarded as an extension of the Biot- 
Savart law for the temporally harmonic case. 

PERFECT CONDUCTOR IN AN EXTERNAL MAGNETIC FIELD 

If a perfectly conducting body is placed in an external 
incident field E 1 , H 1 , then a current and charge distribution K, 
a is formed on the surface S of the perfectly conducting body, 
such that the interior of the body is free of any field (see Fig 
2) . The field produced by the induced sources K, a is referred 
to as the scattered field E s , H s . It can be calculated using 
equations (17a) and (17b) and is superposed with the incident 
field ES H 1 . Since the total field 

E = E 1 + E 8 ( 18 > 
H = H 1 + H S (19) 
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vanishes inside a perfect conductor, the boundary conditions on 
the inside of the conductor surface S~ are given as : 

Etan = 0 for r e s" (20a) 

for r e S~ (20b) 

H t an = 0 

There is a discontinuity in H tan on crossing the boundary 
face, and the following are satisfied for the outside of the 
conductor surface S + : 

for res (21a) 
for r e S + (21b) 

H tan = K x xx 

In this case, the surface current density K is related to 
the surface charge density o through the two-dimensional 
continuity equation 

Vg - K = -iexi (22) 

The expression V s • K denotes the two-dimensional divergence 
of the surface current density K in the boundary face S . 

If there are a plurality of perfectly conducting bodies in 
the incident field, then the boundary face S is composed of 
different subfaces S k , i.e. 

s = Us k (23) 

k 

The scattered field E s , H s then comprises all the subfields 
originating from the individual bodies, and the boundary 
conditions according to equation (20a) , (20b) , (21a) and (21b) 
are satisfied for all the conductor surfaces. 
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The object of a scattering problem then consists, for given 
incident fields E i , H 1 and while satisfying the boundary- 
conditions (20a, b) , (21a, b), to find the current and charge 
distribution K, o on the surfaces of the conductor. 

The term "current distribution" will be used below to denote 
a coupled current and charge distribution, without explicitly 
referring to the charge distribution. It can at any time be 
derived-using -equation- (22 ),•■- and- can therefore be subordinated to 
the term "current distribution" . 

For the field outside the perfectly conducting body, a 
classical uniqueness theorem can be applied: 

The field distribution in a closed volume region G is 

uniquely specified by the sources which it contains and 

by the behavior of E tan or H tan at the boundary of the 

volume region G. 

As regards the uniqueness of the scattering problems in 
question, it is sufficient to satisfy one of the boundary 
conditions in equations (20a, b) or (21a, b) . 

Evaluation of E tan by combining equation (18) with equation 
(20a) gives 

— n x = n x E 1 for r e S (24) 

Since E tan is continuous on crossing the boundary faces, it 
is in this case unimportant whether the relationship is evaluated 
on the outside of the conductor surface S + or on the inside of 
the conductor surface S~ . The scattered field E s can be 
determined using equation (17a) . If, furthermore, the surface 
charge density a is expressed through the surface current density 
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K by using the continuity equation (22) , this gives the integral 
equation for the electric field in integrodif f erential form: 



nx/J 
S 



-ikR 



4x 



R itoe 0 L R 3 R 2J 



da' = n x E^r) 



for r e S 



(25) 



. . Evaluating the tangential magnetic field strength.. by 
combining equation (19) with equation (20b) or equation (21b) 
leads to the relationships 

-n x H s = n x H 1 for r e S~ (26a) 

K-nxH s = n x H 1 for r e S + (26b) . 

Unlike in equation (24) , it is here necessary to take into 
account the side of the face S where the observer point r lies. 
In order to achieve an integral equation independent of this, the 
discontinuous behavior of H tan on crossing the face S must be 
considered in more detail. 

This finally gives the integral equation for the magnetic 
field of a perfectly conducting body with smooth surfaces: 

- K(r) - ax (f ^ 4" + ^ K(*0 x R da' = n x H x (r) 

2 J J 471 Vr 3 rZJ 

. for res (27) 



Since use of the integral equation (27) for magnetic fields 
is, amongst other things, restricted to bodies with smooth 
surfaces and large radii of curvature, the integral equation for 
electric fields will henceforth be used as the basis in the 
present invention. This being the case, instead of resorting to 
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equation (25) , a slightly modified formulation of equation (25) 
will be used: 



n x 



-ikR 



-ikR 



i2S!ftff £ K (r.)da- i— V f T S V s ' .z(«.)dtf 



= n x 



for res 



(28) 



5=J 



Solution of the four" Maxwell equations' whil'e satisfying the 
boundary conditions (2 0a, b) , (21a, b) thus reduces to analyzing 
the integral equation (28) . The solution function sought here is 
the current distribution K(r) on the surfaces of the body in 
question. 



VECTOR BASIS FUNCTIONS 

Taking into account the continuity equation (22) gives the 
following for the current distribution K, o: 



N 

K(r) = £a n v|/ n (r) (29a) 
n=l 

1 N 

= " — S a n v S " VnW 
n=l 

(29b) 

whereby 

\jj denoting a vector basis function, 

V e denoting a two-dimensional del operator in the face S, 

a n denoting a moment (scalar coefficient) , 

N denoting the number of basis functions, 

n denoting a numerical variable. 
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The use of vector basis functions on triangular elements is 
known from "Sadasiva M. Raoi, Donald R. Wilton, Allen W. Glisson 
Electromagnetic Scattering by Surfaces of Arbitrary Shape. IEEE 
Trans. Antennas Propagat . , Vol. 30, No. 3, pp. 409-418, May 1982 
"Ning Yan Zhu and Friedrich M. Landstofer: Application of Curved 
Parametric Triangular and Quadrilateral Edge Elements in the 
Moment Method Solution of the EFIE. IEEE Microwave and Guided 
Wave Letters, Vol ....3,.. ..No- 3 , -pp. .319-321., .. September 199.3., 
presents vector basis functions for parametric elements. It 
includes the idea of using a linear surface current distribution 
to achieve a piecewise constant charge distribution. Since first 
order vector basis functions extend over two neighboring surface 
elements, up to two basis functions interact per cylinder 
element, up to three per triangular element and up to four per 
rectangular element . A rule for the assignment of basis 
functions to surface elements is as follows: 

One basis function and therefore one degree of freedom of 

the resulting system of equations corresponds to each inner 

edge of a discretization grid. 

As a rule, first -order vector basis functions have the 
following properties: 

a) The resulting" surface current density K is continuous. 

b) The resulting surface charge density o is piecewise 
constant . 

c) K_L=0 at the boundary of the definition region of a 
basis function. 

d) Each basis function satisfies the continuity equation 
(22) individually. 
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e) No singular charges occur, as are unavoidable with 
zero-order basis functions. 

f) Each basis function can, at a large distance, be 
regarded as a Hertzian dipole. 

LINE WEIGHTING METHOD 

In the line weighting method, the boundary condition 
E tan -0 is written in the form of a line integral 

J E • dr = 0 (30) 
Cm 

In this case, the curve relates to the m-th element pair, 
consisting of the elements ^ml / S m 2» and extends in a straight 
line from the centroid of the first element to the middle 



of the common edge and then on in a straight line to the centroid 

.c 
m2 



of the second element rfi-, (see Fig. 3) 



The idea with the line weighting method consists in 
converging the act of taking the gradient of the scalar 
potential cp through a line integral into a discrete potential 
different. The starting point is the integral equation for the 
electric field (28) 

(icoA + V«p) tan = E^ an f or r e S (31) 

with the potentials 

-ikR 

A = ^ff 5 K(r') da' 

4* J ] l R (32a) 



S 



e ikR „ ' / x , (32b) 



<p = - — IT — - — V s ' • K(r') da' 



itaeo47t ^ R 
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If both sides of equation (31) are weighted vectorially over 
the M curves then this gives the equations 



/ (icoA + Vcp) • dr = | E 1 • dr mit m=l, 2, 



,M 



(33) 



-m 



-m 



The integration for V<p which is involved here can now be 
replaced by the potential difference between the two end points. 
•If ■ E' 1 'and" A "are "furthermore taken to first approximation as 
constant within a surface element, then for triangular elements 
this gives: 



10) 



_ml. 

2 



(34) 



with m=l , 2 , 



,M, 



c c 

The two local position vectors djni and ^^12 are 
represented in Fig. 3. 

One advantage of the line weighting method consists in 
circumventing the pronounced singularities from equation (25) , 
which significantly simplifies the calculation. 

SETTING UP THE LINEAR EQUATION SYSTEM 

The following relationships are obtained for the potentials 
from equation (32a) and equation (32b) by using equation (29a) 
and equation (2 9b) 

N — ikR 

A(r) = ^ 2>njJ H'n(r') da' ( 3 5a 

n=l Sn 
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<p(r) = - 



N e -ikR 



ioe 0 47l ^ an If 
n = l Sn 



V s • vj) n (r') da' 



(35b) 



with = Sm ^ Sn2 . 

Thus, and using the relationship I n = a ri/ equation (34) can 
be represented as a matrix equation of the form 



Z • I = u 



(36) 



In this case, the matrix Z = [Zmn] denotes an impedance matrix, 
since it relates the current strengths I = [l n ] to the values 
U = [u m ] . Similarly, the inverse matrix Z _1 is referred to as an 
admittance matrix. 

The required moments I n can be determined by direct or 
iterative solution of the equation system (36) . 

After the surface S has been split into the two element 
surfaces S n ]_ and S n 2 , equation (34) gives the following for the 
matrix element Zj^: 



Zmn = ico 



(An + A 21 )^f -(Al2 + A 22 )^ 



+ 912 + <P22 ~ <Pll ~ ^21 



whereby 



J£0 

471 



-ikRq 
// *-z M> n (r')da' 

q R q 

^np 



(37) 



(38a) 



e -ikR q 



<Ppq = - " — ff 

*^ lC08o47t JJ Rg 

Snp 



V s • vj/ n (r') da' 
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(38b) 
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denote the individual contributions to the potential. The matrix 
element Z^, also referred to as a coupling impedance, describes 
the effect of the n-th element pair on the m-th element pair. 
Fig. 4 represents coupling paths between a source element n and 
an observer element m. For example, A 21 denotes the vector 
potential produced by S n2 at the center of S ml . 

One programmed approach for calculating the impedance matrix 
- [Z^] consists- in separate complete evaluation of the individual 
terms according to equation (37) . Since, however, up to 
three basis functions per element are involved in the case of 
triangular elements, another idea would be to calculate the 
impedance matrix Z^ in stages. To this end, the partial results 
according to equation (38a) and equation (38b) are respectively- 
calculated for two triangular elements and, after having been 
provided with suitable factors, are added to the relevant matrix 
elements . 

The right-hand side of the equation system (36) is 
calculated as follows according to equation (34) : 



For an incident plane wave with k as propagation vector and 
E 0 as polarization direction, E 1 can be determined by means of 
the relationship 



U 



m ~ 



p iLc \ d ml _ E i/ X c \ 



(39) . 



E^r) = Eq • e 



,-ik-r 



(40) 
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ELEMENTS WHICH HAVE IMPEDANCE 

If thin layers are modeled using simple surfaces, then it is 
also possible for imperfect conductors with surfaces exhibiting 
resistance, inductance or capacitance to be incorporated quite 
straightforwardly in the existing equation system. For 
conductors of this type, the appropriate boundary condition 
changes from E^ an = 0 to 

Etan = z s • K < 41) 

Because of the skin effect, the surface impedance z s is, 
amongst other things, dependent on a frequency. However, for 
layers whose thickness 1^ is small compared to the skin depth at 
the frequency in question, the current distribution can be 
assumed to be constant over the cross- section. For example, the 
surface impedance z s of a resistive metal layer at low 
frequencies can thus be calculated approximately through the 
relationship 




(42) 



r sp denoting the resistivity of the metal which is used. 

Owing to the new boundary conditions (see equation (41) ) , 
the weighting integral also changes, from equation (30) to 

/ E • dr = Jz s K • dr (43) 
c m c m 

with the quantities occurring on the right-hand side in equation 
(43) to be treated in the form of correction terms in the 
equation system. 
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Fig. 5a shows distributed impedances, involving weighting 
integrals relating to the entire surface element, here 
represented by the surface impedance z s . The correction term 
{i.e. the right-hand side of equation (43)) is 



J z s K • dr = £ l n J z 8 w n ' 3* 



(44) 



C m " Hn 



so that, as a rule, up to nine matrix entries in equation (36) 
need to be corrected per triangular element, and as a rule up to 
16 per rectangular element. 

It is simpler to treat concentrated impedances in the form 
of components, as represented in Fig. 5b. In contrast to 
distributed impedances (see Fig. 5a) , concentrated impedances can 
be represented by a single edge having an impedance. Only a 
single weighting integral and therefore only a single row in the 
equation system (36) are affected. The m-th edge affected by the 
impedance has the following expression 

N 

2 Zrru^n = U m - (45) 
n=l 

the value being defined in equation (37) and the value U m 
being defined in equation (39) . After the additional term has 
been taken over to the left-hand side, this can be treated in the 
form of a corrected diagonal element 

2 mm = Zmm + 2^ (46) 



CONCENTRATED VOLTAGE SOURCES 

In similar fashion to the case of concentrated impedances, a 
single row of the equation system (36) is effected for 



24 



concentrated voltage sources as well: 
N 

ZZmn^ = U m - l4 (47) . 

n = l 

The additional voltage source Um can be taken over unchanged 
to the right-hand side of the equation system (36) , which gives 
the following relationship: 



THE RELEVANCE OF SYMMETRY 

If the arrangement is a symmetric one, then the task of 
solving the equation system can be reduced considerably. In the 
case of a mirror- symmetric arrangement and symmetric excitation 
U m , the number of degrees of freedom can be reduced by half. If, 
however, the excitations are asymmetric, then the excitation 
vector is split into a symmetric part and an asymmetric part . 

The proposed method can also be used to deal with perfectly 
conducting half -planes. 

In the field of EMV simulation, considerations of this type 
can be used, for example, to take account of the metallic floors 
in the measurement rooms, without increasing the number of 
degrees of freedom. 

NUMERIC EVALUATION OF THE INTEGRALS 

In order to calculate the impedance matrix, it is necessary 
to evaluate the integrals from equations (38a) and (38b) 
suitably. If the source point and the observer point are far 
enough away from one another, then it is sufficient to evaluate 
the integrals fully using numerical integration formulae. In 



-«4 



(48) . 
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cases when the source point and observer point are close to one 
another, problems arise, however, from the singularity of the 
integrands, so special analytical methods are used for this. 

The integrals which occur in equations (38a) and (38b) of 
the type 

-ikR 

da' (49a) 



S 

- ikR 

JJ^_ da > (49b) 
S 

have, for R=0 or r=r', a weak singularity which can be split off 
and integrated analytically. The element surface over which the 
integration extends is in this case quite generally referred to 
as S . 

After some rearrangement, the integral from equation (49a) 
becomes : 

.-ikR „ -ikR . 

e — l 



A d£ -R— da ' = Jf< g '- g > e R da ' + 

S S 



-ikR 



( e - e f).fl^__z-i aa . + 



.S 

s 

IG 3 



IG 2 



4-(*-* f ).JJ±da' (50a), 
S 

IG 4 



whereby 
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^ denotes the free corner of the respective element 
(preferably triangle) , 
£ denotes a vector projected into the plane of the triangle 
(indicated by the tilde) . 

A similar rearrangement is found for the integral from 
equation (49b) 

S -S , S , (50b) 

IG 2 IG 4 

The singular integrals IG 3 and IG 4 can be evaluated 
analytically. Numerical integration formulas are used to 
evaluate the integrals IGi and IG 2 , since the integrands are 
continuous and bounded throughout the element surface. 

De 1' Hospital's rule gives the following limit as R— >0 
the integrands from IG 2 : 



e -ikR _ x -j k -e- ikR 
lim = lim = -Ik 



(51) 



The following is correspondingly found for the integrands 
from IGj.: 



lim 
R->0 



(r - s) 



e -ikR _ 1 



= 0 



(52) 



THE CONJUGATE GRAD IENT METHOD 

Further to the possibility of solving the resulting equation 
systems (36) using direct methods, the use of iterative solution 
algorithms is also an option which is favorable in terms of 
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efficiency. The linear equation system to be solved from 



equation (36) is generally represented as 



Ax = b 



rait A e C 



Nxn 



(53) 



C denoting the set of complex numbers. A characteristic of 
iterative solution methods is the generated sequence of 
approximate- solutions x (m) , referred to as iterated -functions . . .If 
this sequence converges with increasing m to the exact solution 
x, then the iteration process can be truncated once the desired 
accuracy has been reached. 

While direct solution methods have an exactly determined 
computing cost, the cost for iterative solution methods depends 
directly on the number of iteration steps needed, and is 
therefore mostly unknown in advance. 

To solve the equation systems, direct solution methods need 
direct access to the individual matrix elements. Conversely, in 
the case of iterative solution methods, the matrix A preferably 
occurs only in the form of matrix/ vector products, and need not 
therefore be given explicitly. Access to corresponding functions 
is sufficient for calculating the matrix/vector products. 
Depending on the iteration method, one or more matrix/vector 
products are calculated per iteration step, so that in the case 
of fully occupied matrices, 0(N 2 ) floating-point operations are 
generally required for this. It is therefore important, when 
iterative methods are used, to keep the number of iteration steps 
required as low as possible. 



The CG (conjugate gradient) method class also includes the 
GMRES method as disclosed by "Youcef Saas and Martin H. Schultz: 
GMRES : A Generalized Minimal Residual Algorithm for Solving 
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by a predeterminable factor e, for example 10~ 4 . If the 
iteration process is started with the null vector as initial 
value, then the iteration can be truncated as soon as the 
inequality 




(57) 



is satisfied. In principle, the truncation criterion should be 
chosen in such a way that the error due to truncating the 
iteration is just less than the discretization error. 

The total number of iteration steps needed to reach the 
truncation criterion is denoted M, with m=l,2,...,M. 

P RECON D ITIONING 

In most problems, the convergence behavior of the iteration 
process can be accelerated by so-called preconditioning. 

Distinction may be made between preconditioning on the left 
and preconditioning on the right. Equation (53) leads to the 
following for the left- transformed equation system 

M _1 Ax = M _1 b <=> Ax = b mit A = M _1 A 

b = irt) (58) 

The matrix M" 1 is referred to as the preconditioner . In the 
case M=A all the eigenvalues of A are, for example, 1 so the 
exact solution is found after just a single iteration. However, 
with preconditioning of this type, the cost of calculating M" 1 is 
just as great as when the equation system is solved directly. In 
principle, it is also true that the effect of the preconditioning 
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is commensurately better if the matrices A and M are similar. 

When using equation (58) , the value B at the start of the 
iteration process must firstly be determined using a simple 
matrix/vector product. On the grounds of efficiency, the matrix 
A is preferably not calculated explicitly. Since, in the CG 
methods, all the matrices occur only in the form of matrix/vector 
products, it is more favorable to evaluate such operations 
sequentially," "that is "to say in' the f orm of "two 'separate 
matrix/vector products. Furthermore, calculating a suitable 
factorization of M will suffice instead of the inverses M" 1 . The 
relevant matrix/vector products must then, however, be replaced 
by suitable back- substitution routines. 

Right preconditioning is carried out in similar fashion to 
left preconditioning: 

A M _1 M x = b <=> A x = b mit A. = AM -1 (59) 

X = M X 

When using equation (59) , the desired solution x must be 
calculated after the actual integration process using an 
additional matrix/vector product according to x = M^x. 

Since the convergence behavior of the iteration process is 
determined by the distribution of the eigenvalues of A, both 
preconditioning variants are equivalent in terms of convergence 
behavior. In comparison with left transformation, right 
preconditioning has, however, the advantage that the residue from 
equation (54) can be adopted unchanged as an error measure. In 
the case of left preconditioning, the residue from equation (54) 
cannot usually be calculated" without additional cost. It is then 
advantageous to use the following form instead of the residue 
established in equation (54) 
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f( m ) = b - AxH = M _1 rH (60) 

In this case, the matrix M" 1 causes distortions which can 
unfavorably affect the truncation criterion. 

JACOBI PRECONDITIONING 

One possible way of preconditioning consists of the diagonal 
matrix 



Mij = { A± 3 i - j 

1° otherwise (61) . 

The use of this so-called Jacobi preconditioning on the left 
or on the right corresponds to a diagonal scaling of A, in which 
the diagonal elements are converted to one by multiplying the 
equation system through by rows or by columns. 

A further advantageous preconditioning method consists in a 
block, variant of the Jacobi preconditioning. In this case, an 
index set of the degrees of freedom 

F = {1,2, ... ,N} 

must firstly be divided according to 

F = U F w 
w 

into pairwise disjoint subsets F w . The number of these subsets 
will be denoted W. The matrix M is then given through the 
relationship 
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for i, j e F w 



Otherwise 



(62) 



In the case of this type of block Jacobi preconditioning, 
the matrices M and M" 1 are block diagonal matrices. The 
matrix/vector product M^x can be calculated efficiently using 
separate factorizations of the individual diagonal blocks. 
• - The convergence behavior is particularly good- if the -index 
sets F w can be assigned to physical subregions . Such subregions, 
in the simplest case cubes, can be produced by stepwise 
subdivision of the body in question in the form of a tree- like 
structure. At each refinement step, the cubes which have already 
been generated are preferably split into up to eight daughter 
cubes with halved side length. Empty regions are usually 
ignored. The coarsest refinement stage, that is to say a single 
block, is denoted stage 0. The stage index increases 
correspondingly by one at each refinement step. 

The degrees of freedom are advantageously permuted before 
the actual iteration process, so that the block structure which 
arises in the equation system can be handled by programming. To 
this end, the indices are rearranged in such away that successive 
degrees of freedom are assigned common subregions. Thereupon, 
the sub-problem assigned to the index sets F w can be extracted 
from A and solved using suitable factorization. The resulting 
factorizations are then combined to form the block diagonal 
matrix M" 1 . 

As an example, Fig. 6 shows the division of the bodywork of 
a car for refinement stages 2, 3 and 4. 

The use of Jacobi preconditioning is not restricted to 
physical subregions of the same size. Especially in the case of 
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discretization with greatly varying element sizes, there is the 
opportunity to subdivide the body into different sizes of 
subregions . Fig. 7 represents adaptive block Jacobi 
preconditioning for the case of a conductor track over an 
infinite earth surface. In order for the current distribution to 
be reconstructed as accurately as possible by the basis functions 
which are employed, the area around the conductor loop is finely 
gridded. It is preferable -to- have -chosen -the subregions in such 
a way that no more than 3 00 degrees of freedom are contained in 
each region. 

THE FFS (FAST FREQUENCY STEPPING) METHOD 

Many practical EMV problems require broadband 
characterization of the object to be examined. Guidelines 
governing both interference immunity and emission extend over 
wide frequency ranges, so that, in order to simulate the 
corresponding measurements, it is necessary to carry out 
extensive series of computations, in which the frequency f , in 
similar fashion to real measurement, is increased from a lowest 
frequency f min for investigation to a highest frequency f^ for 
investigation. So as not to miss any critical frequency range, 
the step size Af which is used should be chosen as small as 
possible . 

One approach for reducing the computation cost while using, 
iterative solution methods results from the observation, that, 
for sufficiently small frequency steps, the solution vectors of 
interest do not differ too much from one another. It is thus 
possible to use the last calculated solution vector as the new 
start value . 
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A further advantageous development consists in determining 
the start value by extrapolation from solution vectors which have 
already been calculated. For the special case of equidistant 
frequencies, the following relationship 

(0) 

xW = Xi_ 3 - 3xi_ 2 + 3xi_ x (63) 

can, for example, be written for quadratic extrapolation of a 
start value x^ 0 ' from the last three solution vectors x^, Xi_ 2 
and x i _ 1 . 

Fast convergence may be observed if one of the equation 
systems is solved using a direct solution method and the 
calculated factorization is subsequently taken as a 
preconditioner for iterative solution of the remaining equation 
system. 

Combining this scheme with the extrapolation approach 
described above gives an effective procedure for good broadband 
analysis of electrodynamic scattering problems. This method is 
referred to below as the FFS method. The algorithm for the 
special case of equidistant frequencies is as follows: 

f • = f • 

• -"-mm 

calculate A and b 

calculate M" 1 : = A (or factorization of A) 
x : = M : b 
if f < f^ 

f :=f + Af 

calculate A and b 

extrapolate a start value x (0) 

solve A M" 1 x = £ iteratively 

x : = M x x (Alg-1) 
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The computation series is, for example, started at the 
lowest frequency f min . After the corresponding equation system 
has been set up, and the matrix A has been factorized suitably, 
the first solution vector is given by back substitution. The 
calculated factorization of A is then taken as a preconditioner 
M 1 for iteratively solving the remaining equation systems. The 
frequency f is for this purpose increased stepwise from f min to 
fmax/ it being necessary "in "the case of "right preconditioning - to 
pay attention that the transformed value x is to be extrapolated 
instead of the actual solution vector x. The actual solution 
vector x can be determined after the iteration process through 
the relationship x: = M~ l x. 

The required computation cost is, in the FFS method, 
determined essentially by the factorization of the matrix M, that 
is to say by the direct solution of the equation system. Using 
the matrix M 1 as a preconditioner, the subsequent iterative 
solution of the remaining equation systems is usually limited to 
a few iterations. 

The sequence in which the individual frequencies are run 
through in the FFS method need not necessarily, as represented in 
the algorithm (Alg-1) , be from f min to f max . . Likewise, variants 
of the method are conceivable which start at the highest 
frequency f max . solve the corresponding equation system directly 
and subsequently reduce the frequency stepwise. 

The preconditioner M" 1 can furthermore be calculated for a 
central frequency f c . The frequencies are then, according to 
choice, run through from f min to f^. from f,,^ to f min or firstly 
from f c to f min and subsequently from f c to f max . The advantage of 
these alternatives consists in the fact that the broadband 
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acceleration effect of the matrix M" 1 can be used in both 
frequency directions. 

MULT I POLE EXPANSION OF THE RETARDED POTENTIALS 

Central to a fast multipole method are methods with which 
the potentials and fields of a given current distribution can be 
characterized approximately by a few scalar coefficients. As 
indicated by the "term " "mult iptfle method", -preferably spherical 
multipole expansions of the potentials or fields are used. 

Since the line weighting method used in this invention (see 
above for details) is based on the retarded potentials A and cp, 
these two quantities are characterized by corresponding multipole 
expansions . 

As an alternative to this, it is possible to make a 
multipole expansion of the electric field strength E. In 
comparison with the multipole expansions of the retarded 
potentials A and <p, the latter does, however, lead to more 
complex analytical relationships which are difficult to handle 
numerically. 

The aim is to find a set of solutions of the homogeneous 
Helmholtz equation 

^V 2 + k 2 j • U = 0 (64) 

in spherical coordinates, 

k denoting the wave number, and 

u denoting a function which satisfies the Helmholtz equation. 

Using these solutions, any function u can be represented in 
the form of an infinite series. The coefficients which occur in 
this case can frequently be determined analytically. 
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If the Helmholtz equation (64) is expressed using the 
spherical coordinates r, S, a represented in Fig. 8, this gives 
the following formulation 

^"all^'ZlJ + ~? — sin + — + JOi = 0 

j?dr^ drJ ^sinsad^ c®) j? Bix ? & fa 2 (65). 

In order to solve this partial differential equation, a 

separation theorem is chosen 

u = U r (r) - U S (3) • U a (a) 

(66) 

in which the factors U r , U a , and U a each depend on only one 
coordinate. After substituting equations (66) in equation (65) 
and multiplying by 

r 2 sin 2 3 
u 

the following equation is found 



sin 2 S 



dr V dr^ Ua d$ \ dS/ U« h« 2 



Ur dr V dr^ Ua dS V dS / U« dot 2 

(67) 

The a-dependent term can be split and, using separation 
constants m, replaced by 

1 d2tJ a _ m 2 

— = -m 

u a da' 2 (68) . 



Equation (67) then gives the following relationship which is 
independent of a 

_i_ _d_ r 2 diM i _d_r sinS dus^ m 2 2 2 _ 

U r dr V dr J U S sin8 dS l 81 * * dS J sin 2 S 

(69) . 
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In order to separate the S- dependency , the corresponding 
terms in equation (69) are expressed using the separation 
constant 1 : 



U$ s 



1 JL( 8 inl> ™f) - -4- « -1(1 + l) 

sin 8 dSV dS ) sin 2 s 



(70) 



This leaves the r-dependent differential equation 



_i_ d_ T 2 dU^ _ 1(1 + 1} + k : 
U r dr v dr / 



2 r 2 = 0 



(71) 



Using the separation theorem from equation (66) , the 
Helmholtz equation can accordingly be converted into three 
ordinary differential equations 



(72a) 



4 f sin s 5H«) 

sin 8 d9 1 dS J 



1(1 + i) - 



m 



sin 2 S 



U$ = 0 



(72b) 



d U a 2 

f- + m^U a = 0 



(72c) 



The equation (72c) is an ordinary differential equation with 
constant coefficients and the solutions 



cc,m ~ e 



ima 



(73) 
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On account of the secondary condition Ua(a+2n) = U a (a), m 
can take only integer values, i.e. meZ (Z being the set of whole 
numbers) . 

Equation (72b) is a special form of the Legendre 
differential equation. Its solutions are the associated Legendre 
functions of the first and second kind 



US.im " ^(cosQ), Q™(cosS) 



(74) 



Almost all Legendre functions have singularities for Q = 0 
and -9 = n, and are therefore preferably ignored. An exception to 
these are the functions pj 11 (cos d) with 1, meN 0 , which are 
related to 1-th order Legendre polynomials P x through the 
relationship 



pf(x) = (-!)-(! - x 2) m ^^_p L(x) 

(75) 

This gives 



pP(cosS) = P^cosS) (76) 

and 

P, m (cos S) = 0 fur m>l 

1 v i (77) 

The zero to second order associated Legendre functions of 
the first kind are as follows: 
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P§(cos S) = 1 P°(cos S) = - cos 2 8 - - 

P°(cos 8) = cos S P^cos S) = -3 sin 8 cos S 

eJ-(cos 8) = - sin 8 p| (cos B) = 3 sin 2 8 

Equation (72a) is a Bessel differential equation with the 
spherical Bessel functions as solution. 



Ur,l ~ jlfrr), yi(kr) 

(78) 



These are related to the known cylindrical Bessel functions 
Jj (x) , Y x (x) through the relationships 

jl(x) = ^ • Jl + 0,5(x) . yi( x ) = ^ • Y l + 0,5( x ). 



(79) 



and can be represented by elementary functions. Thus, for 
example 

. , . sinfkr) /. x cos(kr) 
joM = -5L-2 yo (kr) = — 

jl(kr) . _ cosM (kr) = _ cos(kr) _ sin(kr) 

Dli j (kr) 2 kr ™ ' ( kr) 2 kr 
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Combining the spherical Bessel functions gives the spherical 
Hankel functions of the first and second kind. 

h^(kr) = ji(kr) + i • yi (kr) < 80a > 



4 2) (kr) = ji(kr) - i • Yl (kr) 



(80b) , 

which are equivalent in terms of representing the total solution 
u to the functions j^Ckr), v i (krl Depending on the kind of 
function u, it may be advantageous to favor spherical -Bessel 
functions or spherical Hankel functions. If, for example, the 
function u remains bounded at the origin, then it is sufficient 
merely to apply the functions j]_(kr^ since the functions y^ (kr), 
and therefore the spherical Hankel functions h^(kr) , h^(kr)^ 
become singular for r-*0 . 

The desired solutions U]_ m of the homogeneous Helmholtz 
equation in spherical coordinates are given as 



u lm = bl (kr) • pH(cos *) • e ima (fli) 



b x (kr) denoting a linear combination of 1 st order spherical 
Bessel or Hankel functions, and leNg and meZ. 

For the full solution u, taking account of equation (77) 
gives the series expansion 



oo i 

u = Z 2 c ml ' Una = 
1 = 0 m = -l 



oo i 



~ 2 Z Cna • bi(kr) • pHcos S) 
1 = 0 m = -l 1 ' 



e xma 



with the constant coefficients C]_ m . 



42 



SPHERICAL HARMONICS 



If the angle -dependent terms cos B) and e imot from equation 
(82) are combined, this gives the so-called spherical harmonics 
Y-[ n (S, a) . Given a suitable normalization factor, they are 



defined as follows 



y m (s o) 8 = a + l (1 - ») ! m s) . e ima ( 83 , 

1 v ' \ 47i (l + m) ! 1 

for 0 £ m s 1. Spherical harmonics with m<0 can be defined 
through a symmetry relationship 

Yl - m (S,a) = (-1)^(3, a)* (84) 
The following orthonormality condition is therefore satisfied 

§ y£{*. a) • Yj?'(S, a)* da = S ir • 6^ 

a. (85) 

whereby 

SI denoting the surface of a unit sphere, 
8 denoting a Kronecker symbol defined as 



5 ir = 



1 for 1=1' 

0 otherwise (86) 
The spherical harmonics form on Q a complete set of 
orthonormal functions in terms of the indices 1 and m. One of 
the most important properties of spherical harmonics is the 
consequent fact that any bounded function gO,a) can be 
developed in a series using the Y-j^S, a) spherical harmonics 

9(».«) = Z Z Cml • a) 

1 = 0 m=-l ( 87 ) . 
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Taking into account equation (85), the coefficients C^ m are 
determined through the relationship 



Clm = § g(», a) • Yl m (S, a)* da 
CI 



(88) 



The spherical harmonics of zero to second order are as 
follows : 



y° = -i- 
0 



r-2 



= 1^ 

\327l 



15 _ j _2 q -i2a 
sin <» • e 



— • sin 8 • cos S • e~ ia 
8n 



-ia 



= IX . fl oos 2 , . i) 
V471 V2 2> 



Y? = JA COBS 



= -,/— ■ sin S • cos S • e 



ia 



871 



$ .^-..in..,' 



ia 



V3 271 



15 c -ir, 2 q 

sin <» • e 



ADDITION THEOREM FOR THE GREEN FUNCTION 
The Green function 

e -ikR 

G(r, r>) = with R = ll r ~ r 'll (89) 

describes, to within a constant factor, the retarded potential of 
an oscillating point charge located at a position r' . The 
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development of the Green function in sphe rical solution of the 
Helmholtz equation has the following form: 



e _ I 1=0 m=-l 

-*dk X % 4 2 )(kr') Y™(S', a') j^kr) Y™(S, a)* Vr < V 

1=0 m=-l (90) " 

MULT I POLE EXPANSIONS 

Using equation (90) , the multipole expansion of the retarded 
potential A and <p can be derived. To do this, equation (90) is 
substituted in equations (14a) and (14b) and the order of 
integration and summation is reversed. A rearrangement of this 
type is permissible since the special case r=r' has already been 
excluded in equation (90) . The vector potential A is developed 
by separate application of equation (90) to the three Cartesian 
components A^, Ay und A 2 . 

Different kinds of multipole expansions result depending on 
the position of the observer and source points. The arrangement 
sketched in Fig. 9a, with r>r' , leads to the so-called global 
multipole expansion, in which the external effect of a local 
source distribution G is represented in the form of a multipole 
expansion. The two circles which are sketched with radii d and D 
indicate the boundaries for a near field (less than radius d) and 
a far field greater than radius D) . 

The corresponding series expansions for the retarded 
potentials are given as 





A(r) = -ikno % £ a? m 4 2 >(kr) Yf(S, a)* 



1 = 0 m=-l 
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<P(r) = - — £ Z P? m 4 2 >(kr) <*)* 91a > ' 

80 1 = 0 m=-l 

it being possible to calculate the global multipole coefficients 

q g 
a=[ m and Pfm through the relationships 



*lm = JJJ J ( r ') iO^O *L m (S', «') «V' 
G 

E?m = JIf P( r ') 3l(kr') Y^S', a') dV (91b) 



To simplify the expression, the three scalar multipole 
coefficients of A have been combined to form vectorial 
coefficients a^. 

The case sketched in Fig. 9b, in which the observer points 
are closer to the coordinate origin than the source point, i.e. 
r<r' , leads to a local multipole expansion. In this case, the 
effect of a sufficiently distant source distribution G is 
developed locally. The retarded potentials A and cp become 



A(r) = -ik^ 0 £ J a lm ^l( kr ) Y f( S ' «)* 
1 = 0 m = -l 



cp(r) = - 



ik 



00 

z 

1 = 0 



m=-l 



m 



j x (kr) Yf («. a) 



(92a) 



with the local multipole coefficients 
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G 



(92b) . 




JJj p(r') 4 2) (kr') Y™ (S*, a') dV 



G 

In the scope of a programmed embodiment , the series 
expansions (91a) and (92a) are truncated after a finite number of 
terms. If, for the outer summation index 1=0, 1,....,L, then it 
can be shown by full induction that the resulting L-th order 
multipole expansions always contain (L+l) 2 terms. The total of 
4(L+1) 2 scalar coefficients are therefore required for multipole 
expansion of the retarded potentials A, cp to L-th order. 

According to the value of 1, the most important multipole 
coefficients are customarily referred to as 



According to Fig. 9a and Fig. 9b, the observer region is 
separated from the source region by two spherical surfaces having 
radii d and D respectively. For d<<A, the convergence of the 
multipole expansions is oriented towards the static case and 
depends essentially on the ratio d/D. The series converge better 
the smaller the chosen value of d/D. However, in the case of a 
source region which is larger in electrical terms, i.e. is of the 
order of a wavelength or more, complex interference patterns are 
formed, and L ^ kd is preferably chosen in order to deal with 
them. 



1=0 : Monopole coefficient 



1=1 : Dipole coefficient 



1=2 : Quadripole coefficient 



1=3 : Octopole coefficient 



1=4: Hexadecapole coefficient, etc. 
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MULT I POLE EXPANSIONS OF THE SUBREGIONS 

In the scope of a fast multipole method, the above-described 
subdivision of the body to be investigated is resorted to. 

In each subregion resulting from the subdivision, an 
independent global and local multipole expansion are employed, 
i.e. for each index set of the degrees of freedom P w , there is 
both a series of global multipole coefficients a^ m , p^ m and a 
series of corresponding local multipole coefficients a^ m , Pi m • 
The multipole expansions are in this case made relative to the 
respective cube center r£ . 

If the representation of the current distribution using the 
basis function \j) n (see equations (29a) and (29b) ) is substituted 
in equations (91a,b) and (92a,b) , then the following is found for 
the global multipole expansions in the individual subregions . 

A w (r) = -ikn 0 £ £ ai mfW hJV) if (*, «)* 
1 = 0 m = -l 

*" (r) = " 2 £ P?m,w 4 2 V) !?(*. «)* <«a> 

0 1 = 0 m = -l 

with the global multipole coefficients 

4m, w = S ^ JJ Vn(r') Dl(kf') if («'. a') da' 
neF w Sn 

P?m, w = S *n If V S ' • Vn(r') jl(kr') yf («'. a') da' 

neF w Sn (93b) 

and for the local multipole expansions in the individual 
subregions 
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A w (r) = -ikn 0 2 2 al m w j^kr) Y™(«, a)* 
1 = 0 m = -l 

1=0 m= _i 



with the local multipole coefficients 



(94a) 



a lm, w 



neF w \N W S n 



Pirn, w - 2 *n JJ V 8 ' • v n (r') hfW) Y f ff ') da ' ( Mb) " 
neF w \N w Sn 

The tilde above the spherical coordinates, r, $ and a 

indicates that they are in this case local coordinates relating 

to the respective cube center r^J - When determining the local 

multipole coefficients in equation (94b) , account is taken of 

source regions lying outside the neighboring regions, i.e. 
neF\N w . 

Since the bodies to be investigated in the field of EMV 
simulation (printed circuit boards, wiring, etc.) Preferably lead 
to small subregions, the convergence behavior. of the series 
expansions from equations (93a) and (94a) is determined 
predominantly by the ratio d/D. For the case in which the global 
multipole expansions (93a) are evaluated outside the neighboring 
regions, then the following is satisfied for cubic subregions in 
three spatial dimensions, both for the global and for the local 
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expansions . 

d _ 1 
D >/3 

STABILITY PROBLEMS AT LOW FREQUENCY 

If the body to be investigated is divided up according to 
degrees of freedom or according to basis functions, then 
instabilities in terms of the global multipole expansion of the 
electric scalar potential cp can arise at low frequencies. 

Fig. 10 illustrates the problem with reference to the 
example of a simple conductor loop whose scalar potential is to 
be represented by four global multipole expansions (Fig. 10a) . 
With decreasing frequency, the total charge on the conductor 
surface vanishes since an increasingly uniform current 
distribution is set up. However, on account of the continuity 
equation, the charges in the individual basis functions grow at 
the same time as 1/to (cf equation (29b) ) . So long as the basis 
functions overlap, the charges can cancel one another and meet 
the requirement for a vanishing total charge. However, the 
"stretching" of the basis functions keeps, on the boundary 
elements, charges which can no longer compensate one another 
since they are assigned to different multipoles (see Fig. 10b) . 
If the multipole expansion of cp is truncated after a finite 
number of terms, then this results in truncation errors which 
cause increasing problems as the frequency decreases. 

This situation can be ameliorated by giving the charges the 
opportunity to cancel out before the actual multipole expansions 
are made. One possible programmed embodiment consists in making 
all the global multipole expansions according to elements instead 
of according to basis functions. 
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THE FAST MULT I POLE METHOD (FMM) 

The fast multipole method, hereafter referred to as FMM , is 
disclosed, for example, by V. Rokhlin: Rapid Solution of Integral 
Equations of Classical Potential Theory, Journal of Computational 
Physics, Vol. 60, pp. 187-207, 1985. All ways of implementing 
the FMM for the dynamic case are suitable for scattering problems 
in which subregions are substantially larger than the wavelength. 
In problems of this type, a high order is needed for the 
multipole expansions used. A substantial disadvantage with this 
formalism results from the property that the multipole 
coefficients are not calculated explicitly. Instead of explicit 
calculation, the function to be expanded is "sampled" at discrete 
points on a spherical envelope surface, whereupon the test values 
are directly processed further. Since the multipole terms of 
higher order cannot be suppressed readily in this case, 
interference may arise from aliasing effects which, under certain 
circumstances, can lead to errors in the multipole expansions. 

The method proposed here is suitable, in particular, for 
simulating scattering problems in which the geometrical 
structures of the models and therefore the dimensions of the 
subregions as well are smaller than the wavelength. With these 
assumptions, the interference effects will remain limited, so 
that just a few multipole coefficients will be sufficient for 
describing potentials or field strengths. With the stabilization 
method described here, the method will behave stably even for low 
frequencies . 

It is particularly interesting to use the fast multipole 
method in the field of EMC problems involving radiation. The 
measurements to be simulated are preferably carried out to within 
a few dB . 
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When applying the described iterative solution methods to 
fully occupied equation systems, the matrix/vector product to be 
calculated preferably leads to a complexity of 0 (N 2 ) floating- 



the special structure of the underlying integral equations, using 
the FMM makes it possible to calculate the matrix/vector product 
and therefore also an iteration step in fewer than 0(N 2 ) 
float ing-point operations. 

A condensed representation of the following form will be 
used for the matrix Z 

Z = Z' + LTG (95) , 

Z' denoting . that part of the matrix Z which describes the 

coupling between neighboring subregions 
LTG denoting the remaining part of the matrix Z, which describes 

the coupling between distant subregions, 

i.e. ieF w ,jeF\N w 

In the case when the stabilization measures explained above 
are dispensed with (see section: Stability Problems at Low 
Frequencies) , the Z' matrix is determined according to 



In the stabilized case, additional correction terms are 
involved . 

Fig. 11 represents an observer region BG of the index set 
F W / which is surrounded by neighboring regions NG of the index 
set N w . The neighboring regions NG are therefore directly 
coupled. Distant regions WEG are represented as indirect 



point operations per iteration step. 



In contrast, on account of 




otherwise 



for i e F w , j e N 



w 



(96) . 
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coupling with the observer region BG. Global multipole 

expansions GMPE are made in regions with a cross, while the local 

multipole expansion is made in the observer region BG. 

g g 

The global multipole coefficients aj m w , Pi m ,w ° f the 

g 

individual subregions will be combined below in a vector c • 

1 

Similarly, a vector c contains all the local multipole 
coefficients w , Pi m w • * n the scope of the local multipole 

expansion, the neighboring regions in Fig. 11 are not taken into 
account since the corresponding contributions are already 
contained in the matrix Z' . If the multipole expansions are made 
to L-th order, then the following is found in the case of W 
subregions 

\2 



: 9, c 1 e C 4 *^ 1 ) . 



The vectors c 3 and c 1 are calculated through the 

relationships 

g _ _ (97a) 
C M = G • I 

c l = T . c g = T - G • I (97b) 



i.e. the matrix G makes it possible to determine the global 
multipole expansions in the subregions for a given current 
distribution I. The so-called translation matrix T then 
calculates the local multipole coefficients c 1 therefrom. The 
translation matrix T collects the global multipole expansions 
represented by crosses in Fig. 11 in a local multipole expansion. 
Finally, the local multipole expansions are evaluated at the 
observer points using the matrix L and added to the neighboring 
contributions : 



53 



m 



m 



ZI = Z'l + Lc 1 (98) . 

In the unstabilized case, the matrices G and L are block 
diagonal matrices, the individual diagonal blocks relating the 
multipole expansions of the subregions to the basis functions 
contained therein. Relationships for calculating the matrix 
elements of G follow directly from the global multipole expansion 
according - to equations ■••{ 93a , b) . - T-he- desired -matrix ■ elements 
correspond to the integral terms occurring in equation (93b) . 
Relationships for calculating the matrix elements of L are given 
by substituting equation (94a) into the weighting equation (34). 
In this regard it should preferably be noted that the quantities 
A, cp from equation (34) are to be replaced by the contributions 
of the respective subregions Aw» <p w . 

Elements which have impedance can be taken into account 
according to the procedure described above. Since correction 
terms occur only for neighboring elements, the modifications to 
be made to the matrix Z' remain limited. 

The calculation of the translation matrix T is described 
below. 

TRANSLATION OPERATORS 

Further to the individual multipole coefficients, unique 
characterization of multipole expansion involves specifying the 
coordinate system in relation to which the expansion is made. 
The multipole coefficients are in this case directly dependent 
both on the expansion center (coordinate original) and on the 
definition of the angle values 3 and a (orientation of the 
coordinate axes) . 
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The operators needed in the scope of the fast multipole 
method in order to convert the coefficients of a given multipole 
expansion to a new coordinate system are translation operators, 
since the orientation of the coordinate axes is preserved, in 
contrast to the so-called rotation operators. If the multipole 
expansions are truncated after a finite number of coefficients, 
the required translation operators can be specified in the form 
of- fully -occupied square translation matrices. 

The geometricals mentioned used below are represented for 
illustration in Fig. 12. The position vector of a point P is 
denoted £ in the old coordinate system and r = (r, d, a) 

in the new coordinate system. Furthermore, the old coordinate 
origin 6 is defined uniquely by the position vector 
r' «= (r' , S' , a') in the new coordinate system. The translation 
therefore takes place in the negative r' direction. The 
translation operators considered here are preferably ones which 
convert a given global multipole expansion into a new local 
multipole expansion, that is to say r < r' . 

The starting point is to derive the required translation 
operators from the addition theorem of the Green function (see 
equation (90) ) , which gives the relationship 




1=0 m=-l 



(99) 



for the translation of global monopole terms (l=m=0) . The 



following identities have been taken into account here 




und 




1 
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The fact that equation (99) describes the desired translation can 
be shown in that the left-hand side contains the old monopole 
term which can be replaced by a superposition of new local 
multipole terms. The individual pref actors 



4 2) 0or') Y™(S'. a') 



depend only on the translation vector-r' and can- be calculated in- 
advance and entered in the first column of the translation 
matrix. 

For a programmed embodiment, the required operators are 
preferably determined recursively by successive differentiation 
of equation (99) . 

In order to describe the required differentiation 
relationships, the following differential operators are firstly 
introduced 

a . a 

5 . = — + i • — - 

<5x dy (100a) 

d d 

5_ = — i • — (100b) 

dx. dy 

- d dooc) 

This gives, for example, 

5+ p 0 2 )(kr) Y°(S, a)] = £ k 4 2 )(kr) Y^S, a) ( 101a) 

3_[h ( 0 2) (kr) Y°(S, a)] = k 4 2) (kr) Y^(S, a) doib) 
^[h^kr) Y°(S, a)] = k hf)(kr) Y^S, a) 



(101c) , 
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i.e. first-order multipole terms (dipoles) can be derived 
directly by suitable differentiation of the monopole term. 
Corresponding relationships for multipole terms of higher order 
are derived below. 

The following elementary differentiation relationships can 
be shown for spherical Hankel functions 



a+Jr 1 - 1 " 1 4 2 )(kr)] = kr 1 * 1 h^kr) sin « 
a_[r 1+1 4 2) (kr)] = kr 1+1 4 2 2 1 (kr)sinS 
a«|r* +1 4 2) (kr)J = kr 1+1 hfl^kr) cos S 



-ia 



(102a) 



(102b) 



(102c) 



The corresponding is found for spherical harmonics 



= a™(21 + i) 



*1+1 



-1+2 



(103a) 



a) 



= -b m (21 + 1) 



7 m-l 
'1+1 



r 1 



+2 



(103b) 



Yf (3, a) 
^+1 



- Cl (21 + 1) 



(103c) 



with the constant coefficients 



a l 



= / (l + m + l)(l + m + 2) (1Q4a) 
V (21 + l)(21 + 3) 
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= (1 - m + l)(l - m + 2) (104b) 
1 \ (21 + l)(21 + 3) 

r m _ 1 (1 + m + l)(l - m + l) 

Cl "V (21 + 1)(21 + 3) (104O, 

which result on account of the normalization factor from the 
definition of the spherical harmonics (see equation (83)) . 
'Combining equations (102a-c) with equations (103a-c) gives the 
differentiation relationships for multipole terms of higher order 

a + [4 2 >(kr) Yf (» f a)] = a m (21 + l) y™ + + i(» f a) + 

+ khfl^kr) sin S e ia Y^(S, a) (105a) 

( 2 ) 

a-[hf)(kr) Yf(» f a)] = -b m (21 + i) bJ^l Yj n - 1 1 (S, a) + 

+ khf 2 1 (kr) s in 8 e~ ia yJ^S, a ) ( 1 0 5b ) 

d z 4 2) (kr) Y 1 m (S, a)j = -c m (21 + i) Y 1 m +1 (S, a) + 

+ kh^^kr) cos 8 Y^(B, a) (105c) , 

which can be rearranged further using the recursion relationship 
for spherical Hankel functions 

h^kr) ♦ hfl^r) = (21 + 1) 5^ (106) 
and the recursion relationships for spherical harmonics 
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sin 8 e ict Yf (S, a) = a^" 1 - Y™* 1 ^, a) - aj 1 • Yf^S, a) (107a) 
sin 9 e" ia yf (S, a) = -b^ 1 • Y^S, a) - bf • Y^S, a) (107b) 

cos * Y»(« f a) = cj^ • Y^S, a) - cf • a) 

This finally gives the symmetric dif f erentation relationships 

S + [4 2 >(kr) Y™( S , a)] = a™ k h* 2 ^ Y™* 1 ^, a) + 

+ a"™" 1 k hfl^kr) Yf^S, a) 

a_[4 2 )(kr) Yf(S, a)] = -bj 1 k 4 2 | x Y™ a) - 

- b^* 1 k 4 2 l 1 (kr) Yf.- 1 1 (S, a) 

a z [4 2 \kr) Y» (», a)] = -c™ k 4 2 ii Yf +1 (» f a) + 

+ c l-l kh l-l( kr )Y 1 m _ 1 (S,a) 

which can be used to calculate the desired translation operators. 
To this end, by successive differentiation of equation (99) with 
respect to the coordinates r' , d' , a' , (1+1) 1-th order 
translation operators are inferred from (1-1) 1-th order 
translation operators. It should be noted in this case that 
d = -d'. 



(108a) 



(108b) 



(108c) 
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MULTISTAGE ALGORITHMS 

Further to the described one-stage FMM, there are also 
various multistage variants, in which subregions of different 
size are used. For illustration, Fig. 13 represents relevant 
subregions for the case of a two-stage FMM variant. 

Fig. 13 represents an observer region BG of index set Fw, 
which is surrounded by neighboring regions NG of index set N w . 
The neighboring "regions NG" axe therefore directly coupled. 
Remote regions WEG are represented as indirect couplings with the 
observer region BG. Global multipole expansions GMPE are made in 
regions with a cross, while local multipole expansion is carried 
out in the observer region BG. 

The size of the subregions is set in multistage FMM 
algorithms, preferably proportional to the distance between the 
observer region and the source region. The multipole expansions 
at the various refinement levels can be converted into one 
another using special translation operators, so that a complete 
multipole expansion is advantageously not carried out for each 
subregion of the hierarchic region structure. 

In the dynamic case, multipole expansions of higher order 
are preferably carried for larger subregions. 

Fig. 1 represents a block diagram containing the steps of a 
method for determining an electromagnetic field. In a step la, 
the global multipole expansion is carried out as described above 
for the body which is to be investigated and is subdivided into 
subregions. In a step lb, the local multipole expansion is 
carried out for the body subdivided into subregions. Finally, in 
a step lc, the electromagnetic field of the body is determined by 
superposition from the global and local multipole expansions. 



60 



• * 

The invention is not limited to the particular details of 
the method depicted and other modifications and applications are 
contemplated. Certain other changes may be made in the above 
described method without departing from the true spirit and scope 
of the invention herein involved. It is intended, therefore, 
that the subject matter in the above depiction shall be 
interpreted as illustrative and not in a limiting sense. 
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